#Alexander F. Gazmararian
#afg2@princeton.edu
#January 19, 2024

#Load packages
library(tidyverse)
library(sandwich)
library(lmtest)
library(modelsummary)
library(kableExtra)
library(broom)
library(scales)
library(viridis)
library(here)
library(RankAggreg)
library(forcats)

#Load function to save figures
source(here("code", "fun", "savefig.r"))

#Load data in conjoint format
g <- readRDS(here("data", "qualtrics", "natsurvey_conjoint.rds"))

# Specify model
f.cj <- selected ~ credits_ + target_ + avgprofit_ + addcredit_ + partic_ + projtyp_

# Estimate model
m.select <- lm(update(f.cj, selected ~ .), g)

# Cluster standard errors by respondent
m.cl <- coeftest(m.select, vcov. = vcovCL(m.select, ~ ResponseId))
m.cl <- tidy(m.cl, conf.int = TRUE)
m.cl <- m.cl[-1, ] # Drop intercept
m.cl <- separate_wider_delim(m.cl, cols = "term", delim = "_", names = c("Attribute", "Level"))
# Rename attribute levels
m.cl <-
  m.cl %>%
  mutate(
    Attribute = case_when(
      Attribute == "credits" ~ "Federal Tax Credits for Solar Installation\nBaseline: 0%",
      Attribute == "target" ~ "Target Community\nBaseline: Any",
      Attribute == "avgprofit" ~ "Your Average Solar Yearly Profits\nBaseline: $300",
      Attribute == "addcredit" ~ "Additional Tax Credits for Donation Recipient\nBaseline: None",
      Attribute == "partic" ~ "Community Member Participation\nBaseline: 0%",
      Attribute == "projtyp" ~ "Solar Project Type\nBaseline: Local businesses"
    )
  )
# Create plot
p.select <-
  m.cl %>%
  ggplot(aes(x = estimate, y = Level)) +
  geom_vline(xintercept = 0, lty = "dashed", color = "grey") +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), width = 0, size = 1.25) +
  geom_point(size = 3, shape = 21, fill = "white") +
  geom_text(aes(label = percent(estimate, accuracy = 1, suffix = "")), vjust = -1, size = 3.5) +
  theme_bw(base_size = 12) +
  theme(panel.grid = element_blank()) +
  facet_wrap(~Attribute, scales = "free_y", nrow = 3, ncol = 2) +
  labs(x = "Change in Probability of Support", y = "") +
  scale_y_discrete(labels = label_wrap(28)) +
  scale_x_continuous(limits = c(-.2, .2), labels = percent)
p.select
savefig(p.select, "fig_cj_select", height = 4.75)

# Donation Outcomes
# Estimate model
m.donate <- lm(update(f.cj, donate ~ .), g)
# Cluster standard errors by respondent
m.donate.cl <- coeftest(m.donate, vcov. = vcovCL(m.donate, ~ ResponseId))
m.donate.cl <- tidy(m.donate.cl, conf.int = TRUE)
m.donate.cl <- m.donate.cl[-1, ] # Drop intercept
m.donate.cl <- separate_wider_delim(m.donate.cl, cols = "term", delim = "_", names = c("Attribute", "Level"))
# Rename attribute levels
m.donate <-
  m.donate.cl %>%
  mutate(
    Attribute = case_when(
      Attribute == "credits" ~ "Federal Tax Credits for Solar Installation\nBaseline: 0%",
      Attribute == "target" ~ "Target Community\nBaseline: Any",
      Attribute == "avgprofit" ~ "Your Average Solar Yearly Profits\nBaseline: $300",
      Attribute == "addcredit" ~ "Additional Tax Credits for Donation Recipient\nBaseline: None",
      Attribute == "partic" ~ "Community Member Participation\nBaseline: 0%",
      Attribute == "projtyp" ~ "Solar Project Type\nBaseline: Local businesses"
    )
  )
# Create plot
p.donate <-
  m.donate %>%
  ggplot(aes(x = estimate, y = Level)) +
  geom_vline(xintercept = 0, lty = "dashed", color = "grey") +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), width = 0, linewidth = 1.25) +
  geom_point(size = 3, shape = 21, fill = "white") +
  theme_bw(base_size = 12) +
  theme(panel.grid = element_blank()) +
  facet_wrap(~Attribute, scales = "free_y", nrow = 3, ncol = 2) +
  labs(x = "Change in Donation Allocation", y = "") +
  scale_y_discrete(labels = label_wrap(28)) +
  scale_x_continuous(labels = percent_format(scale = 1), limits = c(-5,5)) +
  geom_text(aes(label = percent(estimate, accuracy = 1, scale = 1)), vjust = -1, size = 3.5)
p.donate
savefig(p.donate, "fig_cj_donate", height = 4.75)

# Treatment Effect Heterogeneity

# Estimate models
m.pid <- list()
m.pid[["Democrat"]] <- lm(update(f.cj, selected ~ . ), g, subset = (pid3 == "Democrat"))
m.pid[["Republican"]] <- lm(update(f.cj, selected ~ . ), g, subset = (pid3 == "Republican"))
m.pid[["Neither"]] <- lm(update(f.cj, selected ~ . ), g, subset = (pid3 == "Neither"))
# Cluster standard errors
pid.cl <- lapply(m.pid, function(x) coeftest(x, vcov. = vcovCL(x, ~ ResponseId))  %>% tidy(., conf.int = TRUE))
# Create dataframe
pid.df <- bind_rows(pid.cl, .id = "Party")
# Drop intercept
pid.df <- subset(pid.df, term != "(Intercept)")
# Create plot
pid.df <- separate_wider_delim(pid.df, cols = "term", delim = "_", names = c("Attribute", "Level"))
# Rename attribute levels
pid.df <-
  pid.df %>%
  mutate(
    Attribute = case_when(
      Attribute == "credits" ~ "Federal Tax Credits for Solar Installation\nBaseline: 0%",
      Attribute == "target" ~ "Target Community\nBaseline: Any",
      Attribute == "avgprofit" ~ "Your Average Solar Yearly Profits\nBaseline: $300",
      Attribute == "addcredit" ~ "Additional Tax Credits for Donation Recipient\nBaseline: None",
      Attribute == "partic" ~ "Community Member Participation\nBaseline: 0%",
      Attribute == "projtyp" ~ "Solar Project Type\nBaseline: Local businesses"
    )
  )
# Create plot
p.select.pid <-
  pid.df %>%
  ggplot(aes(x = estimate, y = Level, color = Party, shape = Party)) +
  geom_vline(xintercept = 0, lty = "dashed", color = "grey") +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), width = 0, position = position_dodge(.8)) +
  geom_point(size = 3, position = position_dodge(.8)) +
  theme_bw(base_size = 11) +
  theme(panel.grid = element_blank(), legend.position = "bottom") +
  facet_wrap(~Attribute, scales = "free_y", nrow = 3, ncol = 2) +
  labs(x = "Change in Probability of Selection", y = "") +
  scale_y_discrete(labels = label_wrap(28)) +
  scale_x_continuous(labels = percent, limits = c(-.3, .3)) +
  scale_color_viridis(discrete = TRUE) +
  geom_text(aes(label = percent(estimate, accuracy = 1, suffix = ""),
                color = Party), vjust = -.75, size = 3.5, position = position_dodge(.8))
p.select.pid
savefig(p.select.pid, "fig_cj_select_pid", height = 6.5)

# Solar Interest
m.interest <- list()
m.interest[["Interested"]] <- lm(update(f.cj, selected ~ . ), g, subset = (buyself_bin == 1))
m.interest[["Not Interested"]] <- lm(update(f.cj, selected ~ . ), g, subset = (buyself_bin == 0))
# Cluster standard errors
m.interest <- lapply(m.interest, function(x) coeftest(x, vcov. = vcovCL(x, ~ ResponseId))  %>% tidy(., conf.int = TRUE))
# Create dataframe
interest.df <- bind_rows(m.interest, .id = "Solar Interest")
# Drop intercept
interest.df <- subset(interest.df, term != "(Intercept)")
# Create plot
interest.df <- separate_wider_delim(interest.df, cols = "term", delim = "_", names = c("Attribute", "Level"))
# Rename attribute levels
interest.df <-
  interest.df %>%
  mutate(
    Attribute = case_when(
      Attribute == "credits" ~ "Federal Tax Credits for Solar Installation\nBaseline: 0%",
      Attribute == "target" ~ "Target Community\nBaseline: Any",
      Attribute == "avgprofit" ~ "Your Average Solar Yearly Profits\nBaseline: $300",
      Attribute == "addcredit" ~ "Additional Tax Credits for Donation Recipient\nBaseline: None",
      Attribute == "partic" ~ "Community Member Participation\nBaseline: 0%",
      Attribute == "projtyp" ~ "Solar Project Type\nBaseline: Local businesses"
    )
  )
# Create plot
p.interest <-
  interest.df %>%
  ggplot(aes(x = estimate, y = Level, color = `Solar Interest`, shape = `Solar Interest`)) +
  geom_vline(xintercept = 0, lty = "dashed", color = "grey") +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), width = 0, position = position_dodge(.8)) +
  geom_point(size = 3, position = position_dodge(.8)) +
  theme_bw(base_size = 11) +
  theme(panel.grid = element_blank(), legend.position = "bottom") +
  facet_wrap(~Attribute, scales = "free_y", nrow = 3, ncol = 2) +
  labs(x = "Change in Probability of Selection", y = "") +
  scale_y_discrete(labels = label_wrap(28)) +
  scale_x_continuous(labels = percent, limits = c(-.25,.25)) +
  scale_color_viridis(discrete = TRUE) +
  geom_text(aes(label = percent(estimate, accuracy = 1, suffix = ""),
                color = `Solar Interest`), vjust = -.75, size = 3.5, position = position_dodge(.8))
p.interest
savefig(p.interest, "fig_cj_select_solarinterest", height = 6.5)

# By conjoint comprehension
m.comp <- list()
m.comp[["High Comprehension"]] <- lm(update(f.cj, selected ~ . ), g, subset = (n_right == 6))
m.comp[["Low Comprehension"]] <- lm(update(f.cj, selected ~ . ), g, subset = (n_right != 6))
# Cluster standard errors
comp.cl <- lapply(m.comp, function(x) coeftest(x, vcov. = vcovCL(x, ~ ResponseId))  %>% tidy(., conf.int = TRUE))
# Create dataframe
comp.df <- bind_rows(comp.cl, .id = "Comprehension")
# Drop intercept
comp.df <- subset(comp.df, term != "(Intercept)")
# Create plot
comp.df <- separate_wider_delim(comp.df, cols = "term", delim = "_", names = c("Attribute", "Level"))
# Rename attribute levels
comp.df <-
  comp.df %>%
  mutate(
    Attribute = case_when(
      Attribute == "credits" ~ "Federal Tax Credits for Solar Installation\nBaseline: 0%",
      Attribute == "target" ~ "Target Community\nBaseline: Any",
      Attribute == "avgprofit" ~ "Your Average Solar Yearly Profits\nBaseline: $300",
      Attribute == "addcredit" ~ "Additional Tax Credits for Donation Recipient\nBaseline: None",
      Attribute == "partic" ~ "Community Member Participation\nBaseline: 0%",
      Attribute == "projtyp" ~ "Solar Project Type\nBaseline: Local businesses"
    )
  )
# Create plot
p.comp <-
  comp.df %>%
  ggplot(aes(x = estimate, y = Level, color = Comprehension, shape = Comprehension)) +
  geom_vline(xintercept = 0, lty = "dashed", color = "grey") +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), width = 0, position = position_dodge(.8)) +
  geom_point(size = 3, position = position_dodge(.8)) +
  theme_bw(base_size = 11) +
  theme(panel.grid = element_blank(), legend.position = "bottom") +
  facet_wrap(~Attribute, scales = "free_y", nrow = 3, ncol = 2) +
  labs(x = "Change in Probability of Selection", y = "") +
  scale_y_discrete(labels = label_wrap(29)) +
  scale_x_continuous(labels = percent) +
  scale_color_viridis(discrete = TRUE) +
  geom_text(aes(label = percent(estimate, accuracy = 1, suffix = ""),
                color = Comprehension), vjust = -1, size = 3.5, position = position_dodge(.8))
p.comp
savefig(p.comp, "fig_cj_select_comp", height = 6.5)

# Create ordered list
ord.list <- subset(g, select = attribute_1:attribute_6)
g_mat <- as.matrix(ord.list)
# Estimate optimal list
optlist <- BruteAggreg(g_mat, 6)
optlist

# Ranking of Preferred Attribute
p.rank <- g %>%
  pivot_longer(cols = c(attribute_1:attribute_6)) %>%
  mutate(name = gsub("attribute_", "", name)) %>%
  mutate(value = factor(value, ordered = TRUE, levels = rev(optlist$top.list))) %>%
  group_by(value, name) %>%
  summarize(n = n()) %>%
  group_by(value) %>%
  mutate(pct = n/sum(n)) %>%
  ggplot(aes(x=value,y=pct,fill=name,label=percent(pct,accuracy=1))) +
  geom_col() +
  coord_flip() +
  theme_bw(base_size = 13) +
  geom_text(size = 4, fontface = "bold", color = "white", position = position_stack(vjust = 0.5)) +
  scale_fill_viridis(discrete = TRUE) +
  scale_y_continuous(expand = c(0,0), labels = percent) +
  scale_x_discrete(labels = label_wrap(20)) +
  theme(panel.grid = element_blank()) +
  labs(fill = "Rank", x = "Percent Ranked", y = "Attribute")
p.rank
savefig(p.rank, "fig_rank_attribute", height = 3)


# Ranking of Preferred Target Community
ord.list <- subset(g, select = targetrank_1:targetrank_5)
g_mat <- as.matrix(ord.list)
# Estimate optimal list
optlist <- BruteAggreg(g_mat, 5)
optlist$top.list

p.community <- g %>%
  pivot_longer(cols = c(targetrank_1:targetrank_5)) %>%
  mutate(name = gsub("targetrank_", "", name)) %>%
  mutate(value = factor(value, ordered = TRUE, levels = rev(optlist$top.list))) %>%
  group_by(value, name) %>%
  summarize(n = n()) %>%
  group_by(value) %>%
  mutate(pct = n/sum(n)) %>%
  ggplot(aes(x=value,y=pct,fill=name,label=percent(pct,accuracy=1))) +
  geom_col() +
  coord_flip() +
  theme_bw(base_size = 13) +
  geom_text(size = 4, fontface = "bold", color = "white", position = position_stack(vjust = 0.5)) +
  scale_fill_viridis(discrete = TRUE) +
  scale_y_continuous(expand = c(0,0), labels = percent) +
  scale_x_discrete(labels = label_wrap(20)) +
  theme(panel.grid = element_blank()) +
  labs(fill = "Rank", x = "Percent Ranked", y = "Target Community")
p.community
savefig(p.community, "fig_rank_community", height = 3)
